Summation by parts and dissipation for domains with excised regions 
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We discuss finite difference techniques for hyperbolic equations in non-trivial domains, as those 
that arise when simulating black hole spacetimes. In particular, we construct dissipative and dif- 
ference operators that satisfy the summation by parts property in domains with excised multiple 
cubic regions. This property can be used to derive semi-discrete energy estimates for the associated 
■ initial-boundary value problem which in turn can be used to prove numerical stability. 

O ' 

o 

(N ■ I. INTRODUCTION 



o 
O 



> 



An important problem in astrophysics is to model in a detailed way the collision of two black holes p], @] ■ This 
requires numerically integrating Einstein's field equations and extracting from the simulations the relevant physical 
information. Unfortunately, it is difficult to obtain numerical solutions of these equations in generic three-dimensional 
settings, especially for long term simulations. Obstacles to this goal are encountered in both the analytical and nu- 
merical arenas. In the analytical one, the formulation of a well posed initial-boundary value problem is not completely 
understood. This includes the definition of proper initial and boundary conditions and the equations determining 
■ the future evolution of the fields. In the numerical arena one seeks to define numerical techniques that allow for 
long term accurate evolutions. This requires the construction of appropriate discrete operators to implement the 
initial-boundary value problem. To date, despite considerable advances in both fronts|23|, the challenge of simulating 
OO ' generic three-dimensional black hole systems remains unattained. 

This article intends to provide some initial steps for setting up numerical techniques suitable to address the numerical 
stability of equations like the ones in question, by extending and devising finite difference techniques to tackle first 
order symmetric hyperbolic problems in non-trivial domains, with numerical stability being guaranteed in the linear 
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case. Furthermore, via a local argument, one can assert that these methods should be useful in evolving smooth 
solutions of quasi-linear symmetric hyperbolic equations as well, as is the case of the full, non-linear Einstein vacuum 
$H . equations when appropriately written 5j. Although the main motivation of this work is to present techniques for 
^P" the simulation of Einstein's equations on domains with excised regions, the techniques here presented are readily 
applicable to any symmetric hyperbolic problem in such domains. Applications of these techniques in a variety of 
scenarios will be presented elsewhere [H H B E3 ■ 

This work is organized in the following way. In sections Illland lTTTI we review some of the issues involved in obtaining 



- r— I 

X 



5_j 

stable numerical schemes through the energy method, to set the stage for the specialized discussions that follow. In 



scction lrvl the main new results of the article are presented: we derive three dimensional difference operators satisfying 
summation by parts for non-trivial domains which enable one to obtain energy estimates. We further introduce 

dissipative operators which do not spoil these estimates. To complete a stability proof, one needs to impose 
boundary conditions without affecting the semi-discrete energy estimates obtained thanks to the use of the constructed 
dissipative and difference operators. One way of doing so is by imposing boundary conditions through an orthogonal 
projection, as done by Olsson [T]| . a technique which can even be applied to non-smooth domains. Conclusions are 
drawn in section^ The appendices include technicalities, discussions about Courant limits and further details about 
the implementation of boundary conditions. 



II. NUMERICAL STABILITY THROUGH THE ENERGY METHOD: AN OVERVIEW 

In order to construct stable finite difference schemes for initial-bound ary value problems (IBVPs) associated with 
partial differential evolution equations we will use the method of lines |l2| . This means that we first discretize the 
spatial derivatives appearing in the partial differential equations so as to obtain a large system of ordinary differential 
equations for the grid functions. This system is usually called semi-discrete system. In the following, we assume that 
the solutions of the partial differential equations satisfy an energy estimate which bounds some norm of the solution 
at some fixed time t in terms of a constant C = C(t) which is independent of the the initial and boundary data times 
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the same norm of the initial data and a bound on the boundary data. As we will discuss in more details below, one 
can derive a similar energy estimate for the semi-discrete system if the partial derivatives and the boundary conditions 
are discretized in an appropriate way. Provided that the constant C involved in this semi-discrete estimate can be 
chosen to be resolution-independent for small enough resolution, this immediately implies numerical stability for the 
semi-discrete problem. An important ingredient in the derivation of such estimates is the summation by parts (SBP) 
property which is the discrete analogue of integration by parts. 

There are many possible discretizations of the partial derivatives and the boundary conditions which are consistent 
with the IBVP and which yield a stable semi-discrete scheme, and in general, the constant C will be larger than 
the corresponding constant C of the continuum problem. By carefully choosing the discretization one can achieve 
optimal semi-discrete energy bounds, in the sense that C = C which means that the norm of the solution to the 
semi-discrete problem satisfies the same estimate as the norm of the solution to the analytic solution. Alternatively 
or complementary, one might want to add artificial dissipation in order to control the high frequency modes of the 
solution which are always poorly represented with finite difference approximations. 

Finally, by discretizing the time derivatives one obtains the fully discrete system which is numerically implemented. 
If the semi-discrete system is stable one can show that the fully discrete system is stable as well, provided an 
appropriate time integrator is chosen. This will be briefly discussed later on. 

To summarize, numerical stability in this approach reduces to 

• Formulate a well posed IBVP for the problem one wants to solve at the continuum. 

• Construct difference operators for the domain of interest that satisfy the SBP property. 

• Optionally or even complementary: 

— Construct dissipative operators that do not spoil the energy estimate. 

— Rearrange the semi-discrete equations to achieve optimal energy estimates. 

• Impose boundary conditions without spoiling the semi-discrete stability. 

• Achieve fully discrete stability by appropriately choosing the time integrator. 

The next subsections give an overview of these points. Section ITTT1 discusses some issues that appear in the IBVP 
case, while section IIVI presents results for the case in which there are inner boundaries and excised regions in the 
computational domain. 



There are many definitions of well posedness though roughly, and without going into too many technical details, 
an IBVP is said to be well posed if 

1. a local in time solution with certain smoothness exists, 

2. the solution is unique, and, 

3. the solution depends continuously on the initial and boundary data of the problem. 

There are different approaches to obtaining well posed formulations of a given problem. A common one is the so 
called energy method [l2| where one seeks an energy norm which has the form 



where u(t, .), the solution of the IBVP at some given time t, lies in some Hilbert space with scalar product (.,.). 
In many physical situations, this Hilbert space can be motivated by the existence of a conserved energy. In general 
however, the expression £ does not need to coincide with any physical energy. For first order symmetric hyperbolic 
linear systems, the Hilbert space can be taken to be the space of square integrable functions over some domain Q, 
with scalar product 



A. Semi-discrete stability 



1. Well posedness 



S(t) = \\u(t,.)\\ 2 = (u,u), 



(1) 




(2) 
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where the symmetrizer H is a symmetric positive definite matrix-valued function on Q which is bounded from above 
and away from zero, and u ■ v denotes the standard scalar product between the vector-valued functions u and v. In 
this case, and for suitable boundary conditions, one has an a priori estimate of the form 

£ (t) < C(t)£ (0) + G(t), t>0, (3) 

where C(t) is independent of the initial and boundary data, and G(t) bounds the boundary data. This means that 
the solution at time t can be bounded by the energy at t = Q and a bound on the energy pumped in through the 
boundary of the domain. In particular, for each fixed t > 0, small variations in the data result in small variations of 
the solution. 

As an example, consider the initial-boundary value problem 

d t u(t, x) = —d x u(t,x), x>0, t > 0, (4) 
u(0,x) = f(x), (5) 
u(t,0) = g(t), (6) 

with smooth initial and boundary data, f{x) and g(x), respectively, that satisfy appropriate compatibility conditions 
at (t, x) — (0, 0). We also assume appropriate fall-off conditions at x = oo. Taking a time derivative of the energy 
with H = 1 and using the evolution equation |0J we obtain 

^-£ = -(u, d x u) - (d x u,u) = u 2 (t, 0), (7) 

where we have used integration by parts in the last step. Using the boundary condition JHJ) and integrating one 
obtains the energy estimate © with C(t) = 1 and G(t) = J Q g(r) 2 dT. 

The energy method is not only a valuable tool in studying the system at the analytical level, but it can also be 
used to produce stable numerical discretizations by considering a discrete analogue. Essentially, this involves defining 
discrete operators with which one reproduces, at the discrete level, the steps taken at the continuum when obtaining 
the energy estimate. In what follows we briefly review the steps involved. 

Before we proceed, we mention that the existence of an a priori energy estimate, where the energy norm has the 
form of an integral over a local density expression, is a sufficient but not necessary condition for well posedness 
[T^ |. A different approach for analyzing stability is the Laplace method which gives necessary conditions for well 
posedness. However, the application of this method to obtain sufficient conditions is rather cumbersome and problem 
dependent. On the other hand, the energy method is considerably simpler and so it is preferred when applicable. 
This is the approach we shall follow in this article. It can be applied to a very large set of physical problems. 



2. Summation by parts 

As we have seen, a tool that is used in the derivation of energy estimates is integration by parts. In one dimension 
(ID) it reads 

rb 

v(x) ■ d x u{x) dx + d x v(x) ■ u(x) dx — v(x) ■ u(x)\ a , 

J a 

where a < b. Let N be a positive integer, Ax = (6 — a) /N and consider the grid defined by Xj = a + jAx, j — 0, N. 
The function u(x) is approximated by a grid function (v,j) = (uq, un)- One of the ingredients needed to obtain 
an energy estimate for the semi-discrete problem is to construct difference operators D (approximating the first order 
partial differential operator d/dx) which satisfy a SBP rule ^3] 



{v,D U )l^+{Dv,u)l^=v r u j \^ , (8) 



with respect to some scalar product 



A' 



(u,v)^ N] — Ax ^2 &ij u i ' V J ■ ( 9 ) 
i,3=0 

Here, the weights E = (try ) must be chosen such that the norm E = = (u, which can be seen as the discrete 
version of that defined in Eq. is positive definite. Achievement of SBP in general requires a careful choice of both 



4 



the scalar product and difference operator; for a given choice of scalar product there might not exist any difference 
operators satisfying the SBP property. In the following, we skip the label £ when it refers to the trivial scalar product 
(Jij = 5ij. Also, when there is no need to specify the domain, we skip the superscript [0, N]. 

Having a ID operator that satisfies SBP with respect to a diagonal scalar product £ = (cry) = [SijUj), one can 
construct a 3D operator by simply applying the ID difference operator in each direction. The resulting 3D operator 
satisfies SBP with respect to the scalar product 

(u,v) s = AxAyAz^o-ijkUijk ■ v ijk , 

ijk 

with coefficients (making a slight abuse of notation) er,y/. = OiOjOk This is not necessarily true if £ is not 

diagonal. For this reason, we will only consider diagonal scalar products in this article. In the absence of boundaries, 
SBP reduces to 

(v, Du)£°°' eo '> + (Dv, u)^ 00 '^ = 0, (10) 

provided suitable fall off conditions are imposed. In this case simple ID operators that satisfy SBP with respect to 
the trivial scalar product are the standard centered difference operators D = Dq] for example, 

(D u)j = " J+ ^ A " 3 " 1 (second order accurate) (11) 

(D u) 3 = -"i+2+8"j+i-8u,,--i+ttj- a (fourth order accurate) (12) 

etc. Defining a 3D difference operator by just applying this one-dimensional one in each direction will satisfy SBP 
with respect to the trivial 3D scalar product. 

As a simple application, consider the initial value problem for the symmetric hyperbolic system, 

d t u(t,x) ^ A l d l u(t,x) + Bu(t,x), ieK 3 ,(>0, (13) 
u(0,x)=f(x), (14) 

where u(t, x) is a vector valued function, the initial data f{x) vanishes for \x\ sufficiently large, and where the matrices 
A 1 are symmetric. Assume for simplicity that the A 1 are constant and that B = 0. At the analytic level, the energy 
norm defined by Eq. (JJJ with H — 1 satisfies £(t) — £(0). If we replace di by the centered differencing operator Di 
in the i'th direction, the discrete energy E{t) = \\u(t)\\ 2 with respect to the trivial discrete scalar product satisfies the 
same estimate: E(t) — E(0), and the scheme is stable by construction. 

3. Dissipation 

Even though one can achieve numerical stability through SBP, it is sometimes convenient to add artificial dissipation 
to the problem. One possible reason for doing so is the presence of high frequency modes which - even if they go 
away with increasing resolution - grow in time at fixed resolution. Although this would not be the case in our present 
discussion, it is well known that the addition of dissipation can aid in stabilizing schemes that would otherwise be 
numerically unstable. An example of this would be the case of system that is stable at the semidiscrete level -for 
example, a symmetric system with discretized with operators satisfying SBP- but becomes unstable because the 
region of absolute stability of the time integrator does not contain a neighborhood of the origin along the imaginary 
axis. 

As an example, consider the initial value problem I|13ll4fl . The standard way 12] to introduce dissipation in the 
problem is to modify the right-hand side (RHS) of the equations 

A { diU + Bu^ A % D lU + Bu + Qf>u + Q ( d v) u + Q d z) u, 

where Di is a differential operator satisfying SBP and where Qd is a differential operator that vanishes in the limit of 
infinite resolution, such that the consistency of the scheme is not altered, and such that 

(u, Qd.u)T, < for all u , (15) 

with some inner product £ for which a discrete energy estimate holds (for example, the one with respect to which SBP 
holds). The dissipative property fTB")! ensures that the discrete energy estimate is not spoiled and might be useful to 
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stabilize an otherwise unstable scheme. Furthermore, Qd can be chosen such that it controls spurious high frequency 
modes in the solution. 

As mentioned above, in the absence of boundaries, standard centered derivatives satisfy SBP with respect to the 
trivial scalar product. Similarly, the Kreiss-Oliger dissipation 

Q d = a{-l) r - 1 (Ax)*- 1 D r + D r _, <J>0, (16) 

where 

D +Uj = , D_u = U -^±± , (17) 

Ax Ax 

denote the one-sided difference operators, satisfies the dissipative property (DP) with respect to that scalar product. 
Furthermore, if the accuracy of the scheme without artificial dissipation is q, choosing 2r — 1 > q does not affect the 
accuracy of the scheme. Notice also that with a slight abuse of notation we have used a to denote the dissipative 
parameter; this should not cause confusion with respect to cr, used in the context of weights of scalar products, since 
the latter have subindices. 



4- Optimal energy bounds 

So far, we have discussed how to obtain energy estimates for the semi-discrete problem by constructing finite 
difference operators that satisfy SBP. These estimates imply numerical stability, by providing a discrete analogue of 
Eq. Although the resulting scheme is stable the discrete estimate in principle might not agree with that found 
at the continuum; that is, the constant C(t) appearing in the discrete analogue of Eq. © might be larger than C{t). 
If one can show that C(t) — C(t) we say that the scheme is strictly stable[2$. The solution to non-strictly stable 
schemes can have unwanted features such as artificial growth in time of the errors. In the limit Aa; — > these features 
would disappear; but one would like to avoid or minimize them even at fixed resolution. 

For first order symmetric hyperbolic linear systems, like the one described by Eqs. (|13I14|) . one can achieve strict 
stability [Tlj by rewriting the partial differential equation in skew-symmetric form: 

d t u = ^AWiU + \d % (A 1 u) +(b- \d % {A 1 )^ u . (18) 

One can show that the resulting scheme is strictly stable with respect to the energy E(t) = \\u(t)\\ 2 defined by the 
trivial discrete scalar product. In contrast to this, the simple discretization dtu — A l DiU + Bu, while yielding a 
stable scheme, does not necessarily yield a strictly stable scheme if the matrices A 1 are not constant. Strict stability 
is particularly useful if the formulation at the continuum admits a sharp estimate. In this case, the construction of 
a strictly stable scheme can be exploited to rule out artificial growth of the solutions @- Consider, for example, a 
system with time-independent coefficients, 

d t u = A i (x)d i u + B(x)u. (19) 

Let H be a symmetrizer, i.e. a symmetric positive definite matrix-valued function H(x) which is bounded from above 
and away from zero and which is such that the matrices HA 1 are symmetric. If H can be chosen such that 

di(HA i ) = HB + (HB) T , (20) 

one can show that the energy given by Eq. is conserved, and that by rewriting the semi-discrete equations as 

d t u = ^A*D lU + ^H~ 1 D i {HA i u) + ^B — ^H~ 1 d l (HA i 

the semi-discrete energy does not grow either. Therefore the numerical solution will not have spurious growth even 
at a fixed resolution. In principle, there is no reason why such a symmetrizer should exist. However, on physical 
grounds, there should be one whenever there is a well defined local energy density. 

The usefulness of having a conserved energy at the semi-discrete level is discussed in detail in Ref. . 
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B. Fully discrete stability 

Proceeding through the steps above described A 1111 A 2|) . and optionally Ijll A 3III A 4|) . one obtains an energy 
estimate for the semi- discrete ■problem which implies numerical stability for this semi-discrete system. 

However, one is of course interested in the stability of the fully discrete problem. One particular simple but powerful 
approach to achieve this goal is to follow the strategy based on the method of lines and employ a time-integrator 
guaranteed to give rise to a stable scheme. A useful feature of this approach is that one can derive conditions for the 
time integrator that are sufficient for fully discrete stability and independent of the details of the spatial discretization 
QHUHiIIjS' Two ra ther straightforward options are given by third and fourth order Runge-Kutta schemes 0,0>EJ- 
In appendix [CI we present an elementary discussion, for a wave equation in a domain without boundaries, that gives 
an idea of possible values of the Courant factor. For more complicated equations and domains see Ref. hj. 

So far our discussion has ignored the presence of boundaries. However these are unavoidable in most problems of 
interest. In the next sections we discuss how to modify the previously described techniques such that one recovers 
an energy estimate in the presence of boundaries. For the sake of clarity in the presentation, we first concentrate on 
simple domains which do not posses inner boundaries and then consider the case of domains with holes. 



III. NUMERICAL STABILITY FOR IBVPS 



We here discuss appropriate finite-difference schemes for simple 3D domains. Consider the following IBVP in the 
cubic domain = [0, 1] x [0, 1] x [0, 1], with maximally dissipative boundary conditions (r denoting its (non-smooth) 
boundary) , 

d t u(t,x) = A i d i u(t,x)+Bu(t,x), xeft,t>0, (21) 
u(0,x) = f{x), xen, (22) 

w+{t,x) = S{x)w-{t,x) + g(t,x), xeT,t>0, (23) 

where the matrices A 1 are symmetric and w±(t, x) represent the in- and outgoing characteristic variables with respect 
to the unit outward normal to the boundary. It is assumed that the coupling matrix S(x) does not depend on time 
and is small enough so as to imply an energy estimate 0] . It is also assumed that the initial and boundary data are 
smooth enough and satisfy compatibility conditions at the intersection between the initial surface and boundaries. 



A. Summation by parts 

If one discretizes the RHS of equation (|21|l according to 

d t u{t) = A i D i u{t,x) + Bu(t,x), xen,t>0, (24) 

where the discrete derivative operator Di satisfies SBP, one will obtain an energy estimate, modulo boundary contri- 
butions that appear after SBP. As already mentioned, a simple strategy to construct 3D operators satisfying SBP is 
to use in each direction a ID operator satisfying SBP with respect to a diagonal scalar product, and this is what we 
do in what follows. 

Existence pro ofs of high order difference operators and scalar products satisfying SBP in ID with boundaries can 
be found in Explicit expressions for these operators can be found in |12l |2C| (some of these operators have 

non-diagonal associated scalar products, and therefore their use beyond one dimensional cases does not guarantee 
numerical stability). In the interior the derivatives are approximated by one of the centered finite difference operators 
mentioned in the previous section, and the operator is modified near the boundary points. The simplest example is 

Duq = D + uq , Duj = DqUj (j = 1, . . . , N — 1), Dun = D-Un , (25) 

where the operators Do, D + , D_ are defined in Eqs. (f 1 If) and (|17|l . The operator (|25J) is second order accurate in 
the interior and first order at the boundaries. It satisfies SBP, Eq. JSJ), with respect to the diagonal scalar product 
Co = o~n = 1/2, dj = 1 for j = 1, . . . , N — 1. For high order difference operators, the operator and the weights in the 
scalar product may needed to be modified at points that are near the boundary points as well. 

If the boundary is not smooth and possesses corners and edges the prescription remains the same, provided the 
vertices and edges are convex. Since in the diagonal case the 3D scalar product is just the product of the ID one, the 
weights aijk are 1x1x1 = 1 in the interior points, and in the outer boundary: 1 x 1 x 1/2 = 1/2 at the boundary 
faces, 1 x l/2x 1/2= 1/4 at the edges and l/2x l/2x 1/2 = 1/8 at the vertices. The presence of concave edges and 
vertices requires modifications, which need to be treated carefully. These will arise, for instance, when considering a 
computational domain with an interior excised region. We discuss this case in section IV. 
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B. Artificial dissipation 



Having introduced the simple difference operator (|25|1 that satisfies SBP in a cubic domain with respect to the ID 
scalar product diag(l/2, 1, 1, . . ., 1, 1, 1/2), we now construct some operators that satisfy the dissipative property with 
respect to the same scalar product. We start with the dissipation operator (|16|l for r = 1 but redefine a and rewrite 
the operator as 

Q d = a(Ax) s D + D- . (26) 

In the absence of boundaries, it satisfies the DP for all s = 1, 2, 3.... Notice that unless s = 1, the dissipation parameter 
a is not dimensionless. As we will see later, this redefinition is convenient when boundaries are present. The goal is 
to define Qd through Eq. i|26|) for i = 1 . . . N — 1, and to extend its definition to i — 0, N such that the DP holds. 
Through a straightforward expansion one gets 

(u, Q d u) [ £ N] = ^-u Q Q d u + a(Ax) s (u, D + D-u) [1 ' N ~ 1] + ^-u N Q d u N 

= -<j(Ax) s (H^ull^) 2 + ^u Q [Q d u - 2a(Ax) s - 1 D + u a ] 
Ax 

+— u n [QdU N + 2<7{Ax) s ~ 1 D^u N ] , 

where the second equality comes from Eq. (|A4|) of appendix lAl 

Certainly, there are many possible definitions of Q d at the boundary points that imply the DP, Eq. I|15l) . The 
simplest one is 

Q d u Q = 2a(Ax) s ' 1 D + u , 

Q d u 3 = o{Ax) s D+D_u 3l (j = 1 ... N - 1), (27) 
Qd^N = -2a(Ax) s ~ 1 D_u N . 

Note that this dissipative term vanishes in the limit Ax — > only if s > 1. On the other hand, in those cases, the 
amplification factor of the high frequency modes depends on the resolution already in the absence of boundaries. 
That is, consistency requires s > 1, and in those cases the amplification factor converges to the non-dissipative one 
as Aa; — * 0. This means that the dissipative operator just constructed cannot be expected to cure difference schemes 
which are unstable. It can only be applied to systems where the amplification factor does not become greater than 
one in the limit of high resolution; that is, to schemes which are already stable in the absence of dissipation. We now 
improve on this. 

In the absence of boundaries, the standard way of correcting this is by considering dissipative operators of the 
form Q16[l with 2r — 1 > q, which satisfy the DP and do not change the accuracy of a scheme which uses g'th order 
accurate difference operators. Therefore, we now look for corrections to the operator at and near the boundary 
points. For simplicity, let us consider the case of a dissipative operator for an otherwise second order scheme; that is, 
assume we are using difference operators of order two in the interior, and first order at boundary points. Then, for 
i = 2 ... N — 2 we define Q d through equation (| 1 with r = 2, that is, we set 

Q d = -a{Ax) s D 2 + D 2 _ , (28) 

where we have redefined tr, as before. The modification at and near the boundary points is motivated by a calculation 
that is similar to the one presented above: 

(u,Q d u)^ = ^Axu Q d u + Ax Ul Q dUl - a(Ax) s (u, D 2 + D 2 _u) [2N2] 

+ Axu N ^iQ d UN-i + -Ax u N Q d u N (29) 



Ax 



2 

\s-l n2 , 



uo ( —QdU + cr(Ax) s D + u 



- ui [AxQdUi + cr(Ax)^ 1 (D\ - 2D+D-) u x ] 
-a(Ax) s (JI.D+.D^III 1 ^- 1 !) 2 

- ujv-i [AxQ d u N -i + a(Ax)^ 1 (D 2 _ - 2D+DJ) u N -i] 
'Ax 



un ( -fQdUN + viAxf- 1 D 2 _u N ) , (30) 
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where in the last equality we have used properties listed in appendix 1X1 

As before, there is more than one way of satisfying the DP, the simplest one being: 

Q d u = -2<7(Ax) s - 2 D 2 + u , 

Q dUl = -a(Ax) s - 2 {D 2 + -2D+D-) Ul , 

Q dUj = —a(Ax) s (D + D-) 2 Uj , j = 2, . . . , N - 2 , 

QdUN-i = -<j(Ax) s - 2 (D 2 _ -2D + D_)u N - 1 , 

QdU N = -2a(Ax) s - 2 D 2 _u N . (31) 

For s — 3 the high frequency modes are now damped in a resolution independent way while the dissipative term 
vanishes in the limit Ax — * 0. In this case, the dissipative operator constitutes a third order correction in the interior 
and a first order one at the boundary points and the two grid points next to it. Adding this operator to the RHS 
of the equations does not affect the consistency of the scheme, but reduces in one the order of accuracy at the grid 
points 1 and N — 1 (recall that the difference scheme was already first order at the boundary points, so the accuracy 
there is not changed). 



C. Boundary conditions 



Finally, the maximally dissipative boundary conditions Q23JI are implemented by multiplying the RHS of H24J) from 
left by an appropriate projection operator |ll| P. For homogeneous boundary conditions, P is the projection on the 
space of grid functions that satisfy the boundary conditions (|23H with g = which is orthogonal with respect to the 
scalar product for which SBP holds. The orthogonality of P makes sure that the energy estimate is not spoiled, 
and hence numerical stability still follows. Furthermore, if P is time- independent, a solution to the resulting semi- 
discrete system automatically satisfies the boundary conditions. Inhomogeneous boundary conditions are discussed 
in appendix IdI 



IV. SBP AND DISSIPATION FOR IBVP IN DOMAINS WITH EXCISED REGIONS 



In the previous section we restricted our discussions to computational domains without inner boundaries. Now we 
extend those results to non-trivial computational domains. To simplify the discussion we restrict ourselves to the 
case where a single cubic box is cut out of the computational domain. However, it is straightforward to see that the 
difference and dissipative operators constructed in this section can also be applied to a domain with multiple cubic 
boxes excised from the computational domain, such that SBP and the dissipative property (DP) are satisfied. 



A. Summation by parts 

We will only analyze the case in which the second order centered difference operator Dq, defined in Eq. (|ll|l . is 
used in the interior. The possibility of constructing higher order accurate operators is discussed in appendix iBl The 
modification of D and its associated scalar product at the outer boundary points has already been discussed in the 
previous section. 

In order to define the scalar product and the difference operators at the excision boundary such that SBP holds, 
we restrict ourselves, without loss of generality, to an arbitrary line which is parallel to the x axis. For points on 
such a line we assume that the difference operator in the x direction depends only on neighboring points on this line 
(as we will see next, this assumption suffices for our purposes). In order to label the grid points on such a line we 
only need the index i and assume that there is an inner boundary at, say, i = 0. We will use centered differences at 
the neighbors i = 1 and i = — 1 (this will also turn out to be sufficient) but leave the scalar product at these points 
undefined for the moment. That is, <7ijk at points i = — 1, 0, 1 takes some values a, (3, 7, respectively. Note that since 
the difference operators at the outer boundaries have already been introduced in section ITTT1 for the sake of clarity 
and without loss of generality, we will now ignore this outer boundary and consider the case where only an inner 
boundary is present. Therefore, in what follows, the subindices range from +00 to —00. 

By splitting the scalar product in the different sub-domains one gets (omitting the sum over j and fc) 

(u, Dv) + (v, Du) = a f(u, Dqv)^ 00 ^ + (v, Dqu)^ 00 ^ + 
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FIG. 1: The figure illustrates the set of grid points belonging to a z = const, plane that passes through the interior of the 
excised cube (left panel), or cuts the top or bottom of it (right panel, also shown in Fig|5|). The weights aijk that appear in 
the scalar product are 1 at interior grid points (circles), 1/2 at boundary faces (squares), 1/4 at convex edges (hexagons), 1/8 
at convex vertices, 3/4 at concave edges (triangles) and 7/8 at concave vertices (pentagons). 



FIG. 2: Examples with interior points at the left and right of an edge, respectively, and to the left and right of a vertex, 
respectively (sweeping from left to right in the figure). 



(3 Ax (u a Dv + v Du ) + 7 ((«, A>«) [1,oo) + («, Aju) [1,do) ) (32) 

ot 7 \ 

= 2^ v ° u - 1 + u o u -i) + P^x {u Dv + v Du ) - -(uqVi + VqUi) . 

We now assume that the modified difference operator at i = has a stencil of width three; that is, we assume that it 
has the form 

qui+rua+su-i 

Du = ^ , (33) 

with q, r, s to be determined later. Taylor-expanding Eq. in Ax one obtains that in order for Duq to be 

consistent with d x the conditions 

q + r + s = 0, q- s = l (34) 

must be satisfied. Inserting Eq. Q33JI and these conditions in Eq. I|32(l one sees that in order for SBP to hold, the 
mixed terms of the form vqu±i must vanish. This yields 

a + 2/3s = , 2/3^-7 = 0. (35) 

Equations (|34|) and l|35|l completely fix, at i = 0, the derivative operator (i.e. the coefficients q, r ,s), and the scalar 
product (i.e., 0): 
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7«i + (a - 7)w - otu-x 

Du ° = 1 — i — Va ' 

(a + 7)Ax 

Notice that, unless a — 7 (which in our case would correspond to i = actually not being part of a boundary, as 
we will now see), the modified difference operator (|3fci|> is only a first order accurate approximation of d x , as is the 
case at the outer boundary points. 

We now explicitly write down the scalar product and difference operator just constructed for the different possible 
inner boundary points: 

• Edges: 

By definition of an edge, one neighbor lies in the interior (and has scalar product weight equal to 1) and the other 
lies at a face (and has scalar product weight equal to 1/2). Therefore, Eq. (|36|l implies that (3 = (l+l/2)/2 = 3/4. 
Thus, the expressions for the derivative operator (depending on the face's side) are, 



interior points are (see Fig. |2l 


(a, A 7) 


Du 


to the left of face 


(1,3/4,1/2) 


ui + uq — 2u-i 
3A.t 


to the right of face 


(1/2,3/4,1) 


2iti — uq — U-\ 
3A.r 



• Vertices: 

Similarly, at a vertex on neighbor lies in the interior (and has scalar product weight equal to 1), and the other lies 
at an edge (and has scalar product weight equal to 3/4), therefore (3 = (1 + 3/4) /2 = 7/8 and the corresponding 
derivative operators are: 



interior points are (see Fig. |5J 


(a, A 7) 


-Duo 


to the left of vertex 


(1,7/8,3/4) 


3iti + uq — 4-u-i 
7Ax 


to the right of vertex 


(3/4, 7/8, 1) 


4«i — Uq — 3u-i 

7Ax 



B. Dissipative operator: second order form 



The dissipative operators computed in section lill Bl need to be modified at the excision boundary. We again assume 
that the grid point under consideration at the inner boundary is at index with <tq = (3, and that at the neighbors, 
cr_i = a, ui = 7. We expand the corresponding scalar product: 

(u, Qdu) 1 ^ 00 ^^ + (u,Q d u)^' 0) + (u, QdM)s ,oo) 
a(u, Q d u)(-°°'-V + 0(u, Q d u)W + 7 („, Q dU f<^ 
aa(A X y(u,D + D-u) < >- 00 --^ 
+ (3 Ax u Q d u + 1 a{Ax) s {u, D+D^u)^^ 

a(|| J D_M|| ( - oo ' 0l ) 2 +7(||Z?_w||[ 1 ' oc )) 2 



-a(Ax) s 



+ u [PQ d Ax + a(Ax) s (aL»_ - 7L> + )] u , 



where the last equality is due to Eq. (IA4I) . 

Therefore one (non-unique, but simple) possibility is to choose: 



] d u = -(Ax) s - 1 (jD+ - afl_)iio 



in which case 



(m, Qd.u) = -<j(Ax) s 



(||^||(-^ ]) 2 + 7 (P-«II [1 ' 00) )' 



Making the expressions for the dissipation explicit: 



• Edges: 



• Vertices: 



interior points are (see Fig. EJ 


(a, 0, 7) 


QdUo 


to the left of face 


(1,3/4,1/2) 


Za(Ax) s - 1 (D + - 2D_)u a 


to the right of face 


(1/2,3/4,1) 


|ct(Ax) s - 1 (2D+ - D-)u 




interior points are (see Fig. [3Jl 


7) 


QdU 


to the left of face 


(1,7/8,3/4) 


Za(Ax) s - 1 {3D + -4L>_)w 


to the right of face 


(3/4,7/8,1) 





C. Dissipative operator: fourth order form 

The modification in this case follows along similar lines: 

{u, Q d u)T. 



(u, Qdu)^ °°' 21 + (u, Qdu%_ l ' 11 + +(u, Q d u) [ £' 0] 



+ (u,Q d u) [ ^' 1] + (u,Qd,u) l E 
= a(u, Q d u)(-°°'-y + a(u, Qdu)^ 1 ^ + /3(u, Q d u) [m 
+ j(u,Q d uf>V + 7 (u,Q«i«) [a ' oo) 



[1,00) 



t(Ax) s (u, (D+D_ fu) 



2 \(o°-2] 



+ aAxu^iQdU^i 



+ [3Ax uoQdUo+^AxuxQdUx - ~fa(Ax) s (u, D+D 2 _u) [2 ' oc) 

-a(Ax) s a (llD+D^- 00 -- 1 ^ 2 + 7 (|| J D + Z?_ii|| 

+ u-i [a(Ax) s ' 1 a(2D+D- - D 2 _) + aAxQ d ] it-i 
+ u [pAxQd + ^Axy-^aDl + 7 r>* )] u 
+ m [~fAxQ dUl + a(Ax) s - 1 7 (L>2 - 2D+D-)] , 



where the last equality is due to Eq. (|A3I) . 
Therefore one possible choice is 



Q d u^ 1 = -a(Ax) s - 2 (D 2 _ - 2£>+L>_)m_i , 
a(Ax) s ~ 2 

Qdu = K — '- (aD^+jDDuo , 

P 

Q d ui = -<7(Ax) s ~ 2 (D 2 + - 2D+D-) Ul . 

The case s = 3 should be the preferred one, since then the amplification factor does not depend on resolution. 
Making the expressions for the dissipation explicit: 

• Edges: 



interior points are (see Fig. |5J 


(a, A 7) 


Qdm 


to the left of face 


(1,3/4,1/2) 


-^a(Ax) s - 1 (2D 2 _ + D 2 + )u 


to the right of face 


(1/2,3/4,1) 


-Za(Axy- 1 (D 2 _+2D 2 + )u 
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• Vertices: 



interior points are (sec Fig. [5Jl 


(<x,P,nr) 


QdU 


to the left of vertex 


(1,7/8,3/4) 


-^{Axy-^ADl + 3D 2 + )u 


to the right of vertex 


(3/4,7/8,1) 


-^(Axy-^SDl + 4D 2 + )u 



V. CONCLUSIONS 



We have presented a set of 3D operators satisfying summation by parts which can be used in non-trivial domains. 
The use of these operators to numerically implement first order symmetric hyperbolic systems guarantee stability 
of the semi-discrete system. Additionally, the use of appropriate time integrators and a consistent treatment of the 
boundary provide a systematic way to achieve a stable implementation. Furthermore, in order to rule out artificial 
growth, we have defined dissipative operators which do not spoil the energy growth estimates. 

The usefulness of these techniques will be highlighted in their application to different problems elsewhere. In 
particular, these applications include the evolution of scalar and electromagnetic fields propagating in black hole 
backgrounds in 3D @ ; scalar fields in 2D with a moving black hole Q ; bubble spacetimes which require evolving the 
5D Einstein equations in the presence of symmetries Q , the construction of schemes capable to deal with the axis of 
symmetry in 2D, axisymmetric scenarios |9j and 3D numerical simulations of single black holes in full nonlinear GR 
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APPENDIX A: BASIC PROPERTIES OF FINITE DIFFERENCE OPERATORS 



In this appendix the definition and some properties of the first order 

3 Ax J Ax 

and second order accurate 

DoU > = 2Ax 

finite difference operators are given. The proofs of the following statements can be found in Ref. f° r instance. 
One can show that with respect to the scalar product and norm 



(u,v) [r ^ =J2ujVjAx, (\H\ lr ' s] ) =(w,u) M 



the following properties hold: 

(«, D+v) lr ' s] = -{D-u, »)I r +i.»+i] +u j v j \ s + 1 



s + l 

r ? 



= — (D+u, w)' r ' s ' — h(D+u, D+v)]^] + UjVj\ 

(u,D-v) [r ' s] = -(£)+ti,i;)[ r _i )a _i] +u j Vj\*_ 1 

= -(D-u, v)[ r>s ] + h(D-u, £>_u) [r . s ] +UjVj\*_ 1 , 

(u,D v) [r ' s] = -(D u,v)[ rta ] + -^{ujV j+1 +Uj + iVj)\'_ t . 
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The following equalities are needed in the derivation of the modified derivative and dissipative operators near the 
boundary 

(u,D u)^ = i(^ +1 )£_ l5 (Al) 
{u,D + D_u)^ = - (\\D_u\\t r > s+1 ^ 2 + u s+1 D_u s+1 -u r _ x D + u r _i , (A2) 
(u,(D + D_) 2 u)^ = {\\D + D_u\\^- 1 ^)\^(2D + D.-D 2 _)u s+1 

-^ D - Us+2 + ^ ( 2D+D - - D V ~ ^if D + u - 2 ■ (A3) 

A sketch of the calculation that leads to the last equality is: 

( U) (D + D-)\)^ = (\\D + D-u\\^y 

+ (ujD_D + D_Uj - D_ Uj D + D_u 3 )\ s r +1 (A4) 



= {\\D + D^u\\^f + ~ {uj^iD\uj^i - UjD 2 _Uj) 



8+1 



\D + D_ut-^Y + (2D + D_ D 2 _) u s+1 
J Ax 

u s+2 n 2 , u r-l /on n r>2 \ "r-2 n 2 

" -Ax~ D - Us+2 + ~Ax~ {2D+D - D+) Ur - 1 -Ax~ D+Ur - 2 ■ 

APPENDIX B: ON THE EXISTENCE OF CERTAIN CLASS OF FOURTH ORDER DIFFERENCE 
OPERATORS SATISFYING SBP WHEN EXCISED REGIONS ARE PRESENT 

Here we analyze the question of whether one can construct high order difference operators satisfying SBP in a 
domain with excised cubic regions. More precisely, we seek difference operators that are fourth order in the interior 
and, as in the rest of the article, we assume that the 3D operator is inherited by a ID one with diagonal scalar 
product, that the operator approximating, say, d x , depends only on points to the left or the right. Unfortunately the 
calculation below shows that such an operator cannot be constructed. However, a difference operator that satisfies 
SBP and is fourth order accurate almost everywhere in the interior could be obtained by means of a different approach 
where one decomposes the domain into cubes |22j . 

1. Fourth order accurate operators for domains without excised regions 

The finite difference operator 

(Qu) 3 =D (l- { -^f D+ D_) Uj = -" j+ 2 + ^i-8 Uj -- 1+Uj -- a (m) 

is a fourth order approximation of d/dx. Strand |2Cj showed that there exists a unique second order accurate 
modification of this operator near the boundaries, given by 

(Qu) = (-48u + 59ui -8u 2 - 3u 3 )/(34Ax) , 

(Qu)i = (ua - uo)/(2Aa:) , 

(Qu) 2 = (Suq - 59ui + 59u 3 - 8u 4 )/(86Ax) , 

(Qu) 3 = (3u - 59u 2 + 64u 4 - 8u 5 )/(98Ax) , 

(Qm)jv-3 = (8wjv-5 - 64ujv-4 + 59wjv-2 - 3u N ) / (98Ax) , 

(Qu) N _2 = (8un-4 — 59mat_ 3 + 59mjv-i — 8u N )/(86Ax) , 

(Qu) N -i = (un - HjV-2)/(2Ax) , 

{Qu) N = (3w A r_ 3 + 8uAr_2-59ujv-i + 48u Ar )/(34A2;) (B2) 

that satisfies SBP with respect to a diagonal scalar product, with 

f 17 59 43 49 1 { . 49 43 59 17 

° l ~ I 48' 48' 48' 48' ' '"'' ' '48'48'48'48 
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Next we discuss the question of whether one can modify this operator near inner boundaries, for the case of excised 
cubic regions, in a way such that SBP holds. 

2. Modification of the fourth order operator near an inner boundary 

Let us assume that a 2D domain has a concave corner at the grid point (0, 0), as is the case of an edge at the inner 
boundary of our domain. The scalar product near this point will have the following structure 



1 


1 


1 


1 


1 


1 


1 


1 


1 


1 


1 


1 


4!) 
48 


4!) 
48 


003 


013 


023 


033 


43 
48 


43 
48 


002 


012 


022 


0"32 


59 
48 


59 
48 


coi 


011 


021 


031 


17 

48 


17 

48 


COO 


010 


020 


030 








17 

48 


59 
48 


43 

48 


49 
48 








17 

48 


59 
48 


43 
48 


49 
48 



(B3) 



1 1 



The weights for < i, j < 3 are unknown. We know that they must be positive and symmetric <7jj = aji. We also 
need to compute a second order accurate difference operator. 

Consider one of the rows with 1 < j < 3. Their contribution to the scalar product (u, Du)h must be zero, as there 
are no boundary terms at the continuum. We have 

3 

(u,Qu) = (Tjlu,!)' 4 '!!)'- 00 '- 11 + AxJ2^ijUi(Du)i + {u,D { ^u) [ ^ +co) , (B4) 

i=0 

where is the fourth order accurate difference operator defined in Eq. (|B1|) . 

The contribution from the first and last term of Eq. i|B4|) can be calculated using the properties listed in appendix 
El In general, one has that 

(u, L> 4 u) [r < sl = ^ (8u j+lUj - u j+2 Uj - . (B5) 



We want to solve 



^(8u ti-i - uiu-i - U0W-2) + Axy^ajjUijDuji 



12 

1=0 

- — (8u 4 m 3 - u 5 u 3 - u 4 u 2 ) = (B6) 

for the coefficients of the modified difference operator (Du)i and the weights of the scalar product cr^ for i = 0, 1, 2, 3. 
The operator (Du)i can be at most a 7 point stencil, i.e. 

/n \ 73"i+3 +72"«i+2 +7lUi+l +7o^i +7-lU,_i +7_ 2 "i-2 +7-3""i- 3 

(Du) t = . 

A larger stencil would give rise to terms in (|B4|) that would not cancel. 
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By inspection of Eq. (|B6|) carefully, one sees that there are many coefficients of the modified difference operators 
that must vanish. This reduces the number of parameters from 7 x 4 = 28 to 6. Eq. (|B6|) becomes 



+a 0j u 



3, (I 3 5, . 

--(a + b )u 3 + a a u 2 + I - - -a + -b J Ui 



1 1 15, , 



a 2 u 4 



a 3 u 5 



1 



1 



- — 2ai W*2 + — - + 2ai u — aiw_i 



- - 2a 2 ) u 3 



-- + 2a 2 ) iti - a 2 it 



1 15 l t . 



:a 3 



+&3M1 - o(«3 + &3)«0 



Unless (Tj = 1 there is no solution; to see this it is sufficient to look at the coefficients of uoU- 2 , Uilt-i, u 2 uo 5 ^3^01 
U3U1, U4U2 and U5U3. 

On the other hand, <jj — 1 is in contradiction with the structure of the scalar product under consideration, see 

El. 



APPENDIX C: COURANT LIMITS 



1. Courant limits 



When discretizing a hyperbolic system of partial differential equations with an explicit scheme the Courant- 
Friedrich-Lewy (CFL) condition has to be satisfied in order to have numerical stability. Below we obtain necessary 
and sufficient conditions for the numerical stability of a 3D wave equation, using a standard von Neumann analysis. 
In particular, we want to determine what is the maximum value of A = Af/Aa; that one can use if artificial dissipation 
is added to the RHS. One can use this information as guide for choosing the Courant factor in more general situations. 

We start by considering the scalar advective equation. 



a. Scalar advective equation 

Consider the advective equation u t = au x , with a a real constant. If the spatial derivative is discretized using a 
second order centered differencing operator and the resulting semi-discrete system is integrated using a Runge-Kutta 
time integrator of order p, the solution at a time step n + 1 can be expressed in terms of the solution at the previous 
time step n as 

< +1 = E^( aA *A))X- (Cf) 

This difference scheme IjClf) is stable if and only if it satisfies the von Neumann condition. The Courant limit is 
a\ < \/3 for the p = 3 case (third order Runge-Kutta) and a\ < 2^/2 for the p = 4 case (fourth order Runge-Kutta) . 
These limits change when artificial dissipation is added to the RHS. The fully discrete system becomes 

I 1 

u l +1 =H~ ( aAtr> o ~ <TAt(Ax) 3 {D_D + ) 2 ) s u n k . (C2) 
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The amplification factor depends on the parameters A = Xa and a — a /a. If, for a given value of tr, we numerically 
compute the value of A beyond which the amplification factor becomes in magnitude greater than one for some 
frequency, we obtain the plot in figure 

Courant limits 

for the advective equation u t = a u t 




FIG. 3: The maximum value of the Courant factor that satisfies the von Neumann condition is plotted as a function of the 
dissipation parameter a. 



b. 3D wave equation 



Consider now the 3D wave equation <p tt = 



written in first order form. As before, we use the second 



order centered difference operator to approximate the spatial derivatives and p-th order Runge-Kutta. The finite 
difference scheme is 



«$ 1 = Ei(A*Q)X* 



s=0 



where u 



Vt,<Px,<Py,(Pz) , 



Q = A^D^ + AMD™ + AWdP - aiAxfiD^D^f 
-a(Ay) 3 (D^D ( ^) 2 - a(Azf (D^ } D^f . 
We assume that Ax = Ay = Az and set A = At/ Ax. In Fourier space the difference scheme becomes 

and the amplification matrix is given by 



s=0 



V 



s 








iRx 





s 





iRy 








s 


iRz 


Rx 


iR y 




s 



where R x — Asin(£), R y = Asin(7y), R z = Asin(£) and S = — 16crA(sin(^/2) 4 + sin(?y/2) 4 + sin(£/2) 4 ). The maximum 
values for the Courant factor are plotted in figure 0] 
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Courant limits 

for the 3D wave equation written in first order form 




, i i i i i i i i i i i i i i i i i i i i i 

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 

o 



FIG. 4: The maximum value of the Courant factor for the stability of the difference scheme approximating the first order 3D 
wave equation as a function of the dissipation parameter a. 



APPENDIX D: BOUNDARY CONDITIONS 



Here we discuss how to discretize inhomogeneous boundary conditions of the type l|23|l . Let P be the orthogonal 
projector onto the space of grid functions that satisfy the homogeneous version of the boundary conditions l|23|l. 
By orthogonal we mean that P is Hermitian with respect to the scalar product for which SBP holds. Eq. I|21|) is 
discrctized in the following way: 

d t v = PADv + Bv + (I - P)d t g , (Dl) 

where I denotes the identity matrix on the space of grid functions. Numerical stability for the resulting semi-discrete 
scheme can be shown under additional assumptions on the boundary data (see [HI). 
The construction of the RHS in i|Dl|) proceeds as follows: 

• Discretize the partial derivatives using a difference operators that satisfies SBP. 

• At each Xijk € T 

— Calculate the outward pointing unit normal n and let A n = A % ni denote the boundary matrix. Note that 
the normal at edges and vertices needs special consideration. In the next subsection it is discussed how 
this normal is to be defined in these cases. 

— Multiply the RHS by Q T , where Q is the (orthogonal) matrix that diagonalizes A n , i.e. Q T A n Q = 
diag(A + , A_, 0), where A± is a diagonal positive (negative) definite matrix. Let (W+, W-, Wn) T = Q T ADv; 

— Apply P. Since it is Hermitian with respect to the scalar product for which SBP holds, this projector 
is non-trivial in the sense that, in general, it overwrites the ingoing and outgoing modes. Consider, for 
simplicity, the case in which there is one ingoing and one outgoing mode (or more generally, the case in 
which the coupling matrix between ingoing and outgoing modes is diagonal). Then, we define 



r (new) _ S (qM) (old)-, 1 

+ _ TTs* { + +w - ) + TT^ ; 



w (ne W ) = _}_( SW {old) +W (old) ) _£!_. 

1 ~r o 1 ~r o 

where S is the coupling that appears in the boundary condition w + — Sw_ + g. Note that only when 
5 = does the outgoing mode remain unchanged, otherwise it is overwritten. In the continuum limit, this 
overwriting is, of course, consistent with the partial differential equation. 
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• Go back to the primitive variables v by multiplying by Q. 

1. Corners, edges and vertices 

It is useful to have stability results that allow for cubic domains, since they are used quite often. The corners in 
2D or edges and vertices in 3D deserve special attention. In particular, the normal to the boundary is not defined 
there. The stability of the whole scheme is affected by the treatment of these points. In Ref. [llj it is discussed how 
to control (i.e. how to achieve numerical stability) the boundary terms that appear after SBP at non-smooth parts 
of the boundary. Here we will summarize some aspects of this treatment. 

For simplicity, let us assume that the domain is a 2D square fl = [0, 1] x [0, 1] and that the partial differential 
equation is a symmetric hyperbolic, constant-coefficients one with no principal terms, 

dtu = A x d x u + A y d y u (D2) 

where A x and A v are constant symmetric matrices. If we define the energy norm to be £ = u T u dx dy, the time 
derivative of the energy norm will be given the boundary terms 

j t £ = J u T A*u\%zl dy + J u T Ay U \ V y Zl dx . (D3) 

In particular, in these boundary terms, there is no contribution due to corners, since they constitute a set of measure 
zero. 

The semi-discrete energy estimate obtained by discretizing the RHS with a difference operator satisfying SBP is 
simply the discrete version of (|D3I) : 

UQ j A x u j) 

- ulA^) . (D4) 

One needs to prescribe boundary conditions to the semi-discrete system at the corners. There, the unit vector n is 
not defined. In pd| it is shown how to define n so that a bound to the semi-discrete energy estimate can be given. 
Essentially, the idea is that by looking at the contribution coming from the corner (x,y) = (0,0) in l|D4|) . which is 
given by 

-ctoUqo ( A V AX + &xA y ) u Q 

(note that the "cross-terms" AyA x and AxA v are not a typo but do result from a non-trivial contribution to the 
discrete energy), we see that we need to control the positive speed characteristic variables in the effective direction 
n = —(Ay, Ax)/ A, where A = ((Ax) 2 + (Ay) 2 ) 1 / 2 . For uniform grids with Ax = Ay this effective unit vector at the 
corners of is 

m(0,0) = (-1,-1)/V2, m(l,0) = (+1,-1)/V2, (D5) 
771(0,1) = (-1,+1)/V2, m(l,l) = (+1,+1)/V2. 

This is equivalent to providing boundary conditions as if the normal was at 45 degrees with respect to the faces of 
the cube. Which data to give might be completely or partially determined by compatibility conditions. As discussed 
before, the way in which the boundary conditions are imposed in pd|. including at corners, is through a non-trivial 
orthogonal projector. 

2. A brief discussion about other boundary conditions at corners 

One might wonder whether or not giving boundary conditions at corners along the normal to one of the two faces 
that define the corner yields an energy estimate. That is, if controlling characteristic fields in one of the directions 
(say, the term u T AyA x u ) implies that the other one is also under control because of some compatibility conditions. 
The general answer is no, but it might work in some particular cases. Let us illustrate this with one example. 



— E = ^ Ayr7,:V v ; .1'>Y., 

3=0 

i=0 
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Consider the 2D symmetric hyperbolic system 

/ 



V 




1 

1 



in the domain x > 0, y > 0. The characteristic fields in the x direction are u (eigenvalue 1) and v (eigenvalue — 1). 
In the y direction they are u + v (eigenvalue 1) and u — v (eigenvalue —1). As boundary conditions in the x = and 
y = boundaries use v = su and u — v = r(u + v), respectively, with |r|, \s\ < 1. 
Defining the energy 



£ 



OO pOQ 



(u 2 + v 2 )dxdy, 



taking a time derivative, using the evolution equations and integrating by parts gives (note that in the discrete 
analogue this is where the "cross terms" comes from) 
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At the discrete level one is essentially left with the discrete version of this expression provided operators that satisfy 
SBP are used. In particular, the corner's contribution is 

-E = .,. + {Ay)^(v 2 - u 2 )\ {0fi) - (Ax)(uv)\ m , 

where the dots represent the contributions from the boundary from points other than the corner. One has to control 
the complete corner term as it will not be canceled or changed by any contribution from other boundary points. 
Suppose Ax = Ay. The projector should set, for example, \{y 2 — w 2 )|(o,o) ~ ( uu )l(o.o) to zero. Suppose one does not 
do this but, instead, controls only the ingoing mode in the x direction (that is, v ) by setting it to zero. This does give 
an energy estimate, since 
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E= ... + A l -{~u 2 )\ {m 



which has the appropriate sign (that is, the corner term can be bounded from above by zero). Suppose on the other 
hand that one does not set it to zero but, instead, couples it to the outgoing mode, v = su, with s = ±1. Then 



-E = . . . - sA(u% , 0) 
and the energy estimate follows or not, depending on the sign of s. 
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